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FOREWORD 


This is a progress report on the research oroject "Numerical Solutions 
of Three-Dimensional Navier-Stokes Equations for Closed-Bluff Bodies." The 
period of performance on this research was Januury 1 through December 31, 
1984. The work is supported by the NASA/Langley Research Center through 
Cooperative Agreement NCC1-68, and monitored by Dr Robert E. Smith, Jr., of 
the Analysis and Computation Division (Computer Science and Application 
Branch), NASA/Langley, MS/125. 
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NUMERICAL SOLUTIONS OF THREE-DIMENSIONAL 
NAV IER-STOKES EQUATIONS FOR CLOSED-BLUFF BODIES 


By 

J. S. Abolhassani 1 and S. N. Tiwari 2 
SUMMARY 

With the present facilities at NASA/Langley Research Center, it is 
economically feasible to compute the three-dimensional flow about a complex 
configuration such as closed-bluff bodies (e.g., circular and elliptical 
cylinders) on a flate plate. If the body is sufficiently bluffed, there 
will exist a three-dimensional region close to and about the junction of the 
body with the surface. 

In the present study, the Navier-Stokes equations are solved numerical- 
ly. These equations are unsteady, compressible, viscous, and three-dimen- 
sional without neglecting any terms. The time dependency of the governing 
equations allows the solution to progress naturally for an arbitrary initial 
guess to an asymptotic steady state, if one exists. The equations are 
transformed from physical coordinates to the computational coordinates, 
allowing the solution of the governing equations in a rectangular parallel- 
epiped domain. The equations are solved by the MacCormack time-split tech- 
nique which is vectorized and programmed To run on the CDC VPS 32 computer. 

The todes are written in 32-bit (half word) FORTRAN, which provides an ap- 
proximate factor of two decreasing in computational time and doubles the 
memory size compared to the 64-bit word size. 


Graduate Research Assistant, Department of Mechanical Engineering and 
Mechanics, Old Dominion University, Norfolk, Virginia 23508. 

2 Eminent Professor, Department of Mechanical Engineering and Mechanics, Old 
Dominion University, Norfolk, Virginia 23508. 


1. INTRODUCTION 


The comprehension and analysis of three-dimensional fluid behavior 
around bluff bodies is of considerable importance for flight applications. 
The primary motivation to solve such problems is to understand the basic 
phenomenon of separated flows and to determine the associated forces on the 
object. The static and dynamic features of the unsteady flow field around 
the bluff -bodies are of interest to the designers as well. These flows are 
three-dimensional in the most realistic cases, and can be described by a set 
of Navier-Stokes equations. These equations can be solved nunerically if 
proper grid distributions can be generated. This step (grid generation) is 
the essential first step in solving Navier-Stokes equations for complex 
configurations. Tiiere are presently only a few three-dimensional solutions 
of the Navier-Stokes equations documented, and the obvious reason for this 
is the complexity of the grid arrangements, numerical coding, and lack of 
sufficient numerical resolution (i.e., memory) to generate credible results. 
Most of the available solutions are for simple three-dimensional flow (e.g., 
three-dimensional corner, spherical dome, conical bodies, and spike-nosed 
bodies) which require simple grid arrangements. For many complex geome- 
tries, the constraints imposed by digital computers are less restrictive for 
two-dimensional computation than for three-dimensional ones. One major 
difference between numerical simulation of three-dimensional and two-dimen- 
sional flow is the overwhelming discrepancy in the data base. For a mesh 
with 100 points on a side, a two-dimensional problem requires a memory of 
about 80,000 words (counting 8 dependent variables at each mesh point in a 
typical solution of a compressible flow with arbitrary geometry), while a 
three-dimensional problem requires a memory of about 14,000,000 (14 vari- 
ables) words. The former can be processed easily by any modern scientific 
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computer, however, this is certainly not possible for the la*-.dr case. 

In order to study three-dimensional flows, the Navier-Stokes equations 
are being solved numerically. These equations are unsteady, compressible, 
viscous, and three-dimensional without neglecting any terms. The time de- 
pendency of governing equations allows the solution to progress naturally 
from an arbitrary initial guess to an asymptotic steady state, if one 
exists. The equations are transformed from physical coordinates to the 
computational coordinates, allowing the solution of the governing equations 
in a rectangular parallelepiped domain. The equations are solved by the 
MacCormack time-split technique which is vectorized and programmed to run on 
the CDC VPS 32 computer. The codes were written in 32-bit (half-word) 
FORTRAN, which provides an approximate factor of two decrease in computer 
time and doubles the memory compared to the 64-bit word size. Solutions 
will be obtained for an infinite right circular cylinder on a flat plate at 
high Mach and Reynolds numbers. In this study, the computational planes are 
perpendicular to the axis of the cylinder rather than the flow direction. 
Therefore, updating the boundary and initial conditions become mo^e complex. 

The algebraic method is used to generate the grids and they are dis- 
tributed exponentially in the vertical direction (parallel to the axis of 
the cylinder). Grids are concentrated near solid boundaries and in the 
vicinity of the cylinder-plate junction. The circular cylinder geometry 
would require about 200,000-300,000 grid points and this will occupy 3-4.5 
million 32-bit words of primary memory. It would take about 4-7 seconds to 
complete a time step of explicit MacCormack method. Also, it would take 3- 
4.5 seconds to transfer the restart file from primary memory to secondary 
memory. 
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2. GOVERNING EQUATIONS 


The governiny equations for a thermal fluid system are the conservation 
of mass, momentum, and energy. These equations are developed for an arbi- 
trary region assuming the system is in continuum. Equations of motion for 
viscous, compressible, unsteady, heat conducting flow can be written as: 

Continuity — + V • (pU) = 0, (2.1) 

3c 

Momentum: + v * ( p ^u - T ) = 0, (2.2) 

Energy: — — - + V • (Eu + Q - u • x) = 0. (2.3) 

3t 

V 2 

where E is the total enargy per unit volume given by E = p (e + — + 

<- 

potential energy + ...) and e is the internal energy per unit volume. 
For Newtonian fluid, stress tensor can be related to the pressure and ve- 
locity components as: 


T 


ij 


3U, 3U, 

P « H + u [( _ L + _i ) 

3 Xj 3Xi 


i 

3 ’- 1 3X k 




(2.4) 


This equation is valid under negligible bulk viscosity. For an isotropic 
system, the heat flux in Eq. (2.3) can be expressed in terms of temperature 
gradient (Fourier's law of heat conduction) as: 

q = - K V T (2.5) 

where K is a coefficient of thermal conductivity. A common approximation 

j 
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used for viscosity is based on the kinetic theory of gases using an ideal 
i zed intermolecul ar-forces potential, the relation is: 


where 




T + S 

o 


5 = 198.6° R 
o 

u r = 0.1716 np 


( 2 . 6 ) 


coefficients of thermal conductivity K can be determined from Prandtl 
number 


Pr 

where C y is the specific heat at constant volume and y is the ratio of 
specific heats. 

It is necessary to have a supplementary relation to close the system of 
equations (Eq. (2.1) - (2.3)). By neglecting i ntermolecul ar forces (a ther- 
mally perfect system), thermodynamic properties can be related as: 

P = pRT (2.8) 

where R is the gas constant. Thermally perfect gas assumption allows 
expression of the internal energy (e) as a function of T only [e 1 e(T)]. 
In addition, assumption of calorically perfect gas [e ( 0 ) = 0] a'lows the 
following relation: 
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e = C T (2.9) 

v 

A combination of Eqs. (2.8) and (2.9) results in 

P = pe (ir-1) (2.10) 


These equations (Eqs. (2.1) - (2.3)) are in conservative form. For sim- 
plicity, these equations can be expressed into a compact vector form as: 


where 


U = 


3U 

♦ it 

♦ it 

3H 

+ 

at 

3x 

3y 

3z 


( 2 . 11 ) 


_ 


■ 

P 


PU 

PU 


puu - T + P 



XX 

nv 

, F = 

PUV - T 



xy 

P w 


PUW - T 



xz 

E 


Eu -» q - A + Pu 

A A 





PV 


p w 

PUV - T 


PUW - T 

xy 


XZ 

PVV - T + P 

, H = 

PVW - T 

yy 


yz 

PVW - T 


P WW - T + P 

XZ 


zz 

Ev + - 4>y + PV 


Ew + q z - <(> z + Pw 


For the sake of generality, we can transform these equations from a physical 
domain to a computational domain as 
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( 2 . 12 ) 



3G 3H\ 
3n 3n j 


The transformation coefficients can be computed from a functional re- 
lation between the computational coordinates and the physical coordinates. 


X = x(5,n,5), 

y = y(5,n,5), 

and z = z(5 

,n,5). 

(2.13) 

5 - 5(x,y,z), 

n = n(x,y,z), 

and 5 = 5 ( x 

,y,z) ■ 

(2.14) 


If Eq. (2.14) is known, the transformation coefficients can be computed by 
direct differentiation. If the former relation is not known, after some 
algebraic manipulation, the transformation coefficients can be computed by: 


» ■ 


35_ 

n 

IL 

3 x 

3y 

3 Z 

3n 

3n 

3n 

3 x 

3y 

3 Z 



H 

3 x 

3y 

3 Z 


3x 

9 x 

3 x 

35 

3n 

35 

ay 

Iz 

9y 

35 

3n 

35 

3_Z 

1± 

3_Z 

35 

3n 

35 


(2.15) 
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where 



(i2 

3 Z 
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3x_ 
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35 

3n 

3 C 

3_y 

ll 
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35 

3n 

3 C 

3_Z 

3 _Z 

3_Z 

35 

3n 

3 C 


3x ^3y 3 z 
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3 x ^3y 3 z 
3 c 35 3n 


3y 3 z j 3x ^3y 3z 3y 3z^ 

3C 3n 3n 3; 3n 35 3c 

3y 3z j 
3n 35 


In the present case, the planes of grid are parallel to x-direction thus 
allowing us to write: 


x = x(5 ) 

(2.17) 

y = y(n.c) 

(2.18) 

z = z(n ,c ) 

(2.19) 


This reduces the metric coefficient from n^ne to five non-zero elements. 

3. METHOD OF SOLUTION 

A time marching method is used to compute the solution. This allows us 
to capture the possible transient feature. This method is an explicit 
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second-order accurate time-split predictor-corrector algorithm [3]. The 
governing equations (Eq. (2.12)) are discretized in computational direc- 
tions. In a compact form, they can be expressed as 


u n+1 . „ 

' .J.k 


( L n'“n>] [w] [ L C (4t £>] [W] [ L n (4t n>] 


u. . . 

1 ,J,k 


(3.1) 


where 


At = At 

n ; 



and L r , L , and L are the operators in C, n, snd c directions, 

C n C 

respectively. A time step is completed in this algorithm with the appli- 
cation of each operator applied synmetrical ly about the middle operator. 
For example, operator can be defined as 

L ? (*t ? ) . ^ (3.2) 

where 


Predictor step: 




U i,j,k ' U i,j,k 


At r r 

— ( p i 
AC L 


F, 


i-i) 3- 1 + ( G i * G i-i> r. 1 

3x o y 


+ 



H 




j,k 
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Corrector step: 


u ou ! , = I (u!°. / + 

i » J » ^ 2 \ 1 >J > k 


i,j .k 


At, 


AC 


[■ 




i+1 


- F , ) — i 
' ax 


* < G m - G i> * 1 + < H i*i 



This method has a time' step stability limit, but there is no rigorous sta- 
bility analysis available for this. 'A conservative time step that is com- 
monly used is 


4t < min M + M + M.t c /X7T7 X 

Ax Ay Az \ Ax 2 Ay 2 Az 2 

where c is the local speed of sound. 

In the supersonic region, there exists a large gradient which requires 
a very fine mesh to resolve it. If they are not »esolved, they produce a 
large oscillation which eventually b'lows up the solution. These oscil- 

lations of "low frequency" can be suppressed by adding a fourth order damp- 
ening. A common dampening used is the pressure dampening. This can be 

expressed in physical coordinates as 


-l 

(3.3) 


-°t At 



3 

l v £ ! + c , 

3 2 P 

3U 

3« t 

4P 

3 2 <5 

* 

36 *. 


* = 1,2,3 


(3.4) 


where 6 L = c, 6 2 = n, and 6 3 = c • 
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4. INITIAL AND BOUNDARY CONDITIONS 


In computational fluid dynamics the initial conditions usually corre- 
spond to a real initial situation for a transient problem, or a rough guess 
for a steady state problem. In practice, initial conditions a^e cbtained 
from experiments, empirical relations, approximation theories, or previous 
computational results. An inappropriate initial guess may result in gener- 
ating unphysical ly strong transient waves which propagate through the compu- 
tational region dominating the flow field and eventually lead to a solution 
failure. In general, there are two important requirements that should be 
consiaered in the choice of initial conditions. First, they should be com- 
patible with the fixed upstream boundary conditions. Secondly, the initial 
conditions should be as physically close as possible to the actual nature of 
the How field in the region under study. The former will minimize the 
nunber of iterations required for convergence. An attractive approach is to 
initialize the entire flow field (including the upstream boundary and the 
body surface) with a crude and simple guess (e.g., free stream condition). 
Then, during the course of the computation, both body and upstream boundary 
conditions are changed in a gradual manner to their final values over a 
prescribed number of iterations. The former approach is applied in only one 
step which is equivalent to impulsive initial conditions. 

It is equally important to implement a realistic, accurate, and stable 
method to determine boundary conditions. The application of certain condi- 
tions may cause numerical instability even though the flow is physically 
stable. There are neither mathematical nor physical justifications to im- 
plement a realistic boundary condition. Most of the boundary conditions 
currently implemented are drawn mainly upon intuition, wind tunnel experi- 
ence, and computational experimentation. There are three general types of 
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boundary conditions. They are Dirichlet conditions (specified function 
value), Neumann conditions (specif iced normal gradient), and Robin condi- 
tions (a combination of both). Four important factors should be considered 
in the selection of bundary conditions. They are onvergence, stability, 
computer time, and above all the physical justification. 

For this problem there are five different boundary conditions. They 
are upstream, downstream, lateral, top, and solid boundary. The upstream 
boundary conditions are the undisturbed free stream conoUions and are lo- 
cated at a grid space away from the leading edge, i.e., 


U 



upstream 


(4.1) 


A zero gradient in y-direction (parallel to the primary direction of flow) 
is assumed for the downstream boundary, i.e., 


= 0 (4.2) 
downstream 

The lateral boundaries are located far enough to avoid any influence on the 
interaction region. A boundary-layer profile can be prescribed on the lat- 
eral boundaries. These profiles can be obtained from their corresponding 
points of a flow over a flat pate. Presently, a zero gradient in z-direc- 
tion is assumed for these boundaries, i.e., 


a U 

3y 


15 

az 


= o 

lateral 


(4.3) 
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The cylinder is assumed infinite in height, therefore, the flow at the top 
of the cylinder would be two-dimensional. Consequently, a zero-gradient 
boundary condition is imposed in x-direction for the top boundaries, i.e.. 


ax 


= o 

Top 


(4.4) 


The wall is assumed impermeable and no slip boundary conditions are 
applied, therefore, all velocity components are assumed to be zero. The 
wall is also assumed to have a constant temperature Tw. A zero normal 
pressure gradient is assumed for the solid surface, i.e.. 


3P 

an 


= 0 
sol id 


(4.5) 


This evaluation may appear to be based on the boundary-layer approximation 
(zero normal pressure gradient). In fact, it is a much milder approxima- 
tion, since constant pressure is not applied through the boundary layer but 
over one grid line in the boundary layer. This approximation has yielaed 
stable computations for both non-separated and separated boundary layers 
[ 2 ]. 


5. POST PROCESSING OF DATA AND DATA DISPLAY 
The solution of a three-dimensional flow can produce upto 1,000,000 
data. This vast amount of data needs to be analyzed and necessary informa- 
tion extracted. Storage and manipulation of the huge amount of data becomes 
a serious problem. Color computer graphics provide adequate solutions to 
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the problem. Data which sometimes is over 1,000 pages of computer print-out 
can be compacted and displayed as a digital image on color monitors. Color 
computer graphics offer several advantages for technical presentations, 
documentations, and easy and rapid analysis. A digital image of a flow 
field variable (e.g. density) can be created by transforming node points 
from the computational grids into the object space to the image plane and 
fill the void between the image node point with a color that varies accord- 
ing to the magnitude of the variables. 

6. PRESENT ACCOMPLISHMENTS 

The equations of motion are transformed from the physical coordinates 
to the computational coordinates allowing the solution of the governing 
equations in a rectangular parallelepiped domain. The equations are being 
solved by the MacCormack time-split technique which is vectorized and pro- 
grammed to run on the CDC VPS 32 computer. The codes are written in 32-bit 
(halfword) FORTRAN, which has provided an approximate factor of two decrease 
in computer time and doubles the memory size compared to the 64-bit word 
size. 

The algebraic method is used to generate a 0-type grid around the cyl- 
inder under study, and they are distributed exponentially in the vertical 
direction (parallel to the axis of the cylinder). Grids are concentrated 
near solid boundaries and in the vicinity of the cylinder plate junction. 
Further accomplishments will be reported in a subsequent progress report. 
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